home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / fractal / kaos.lha / eigenlib / hqr2.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-12-15  |  13.9 KB  |  413 lines
  1. hqr2(nm,n,low,igh,h,wr,wi,z)
  2.  
  3. int n,nm,*igh,*low;
  4. double **h,*wr,*wi,**z;
  5. {
  6. int i,j,k,l,m,en,ii,jj,ll,mm,na,nn,itn,its,mp2,enm2,ierr,null;
  7. double p,q,r,s,t,w,x,y,ra,sa,vi,vr,zz,norm,tst1,tst2;
  8. double sqrt(),fabs(),copysign(),cdiv();
  9. char notlas;
  10.  
  11. /*    this subroutine is a translation of the algol procedure hqr2,
  12.       num. math. 16, 181-204(1970) by peters and wilkinson.
  13.       handbook for auto. comp., vol.ii-linear algebra, 372-395(1971).
  14.  
  15.       this subroutine finds the eigenvalues and eigenvectors
  16.       of a real upper hessenberg matrix by the qr method.  the
  17.       eigenvectors of a real general matrix can also be found
  18.       if  elmhes  and  eltran  or  orthes  and  ortran  have
  19.       been used to reduce this general matrix to hessenberg form
  20.       and to accumulate the similarity transformations.
  21.  
  22.       on input
  23.  
  24.          nm must be set to the row dimension of two-dimensional
  25.            array parameters as declared in the calling program
  26.            dimension statement.
  27.  
  28.          n is the order of the matrix.
  29.  
  30.          low and igh are integers determined by the balancing
  31.            subroutine  balanc.  if  balanc  has not been used,
  32.            set low=1, igh=n.
  33.  
  34.          h contains the upper hessenberg matrix.
  35.  
  36.          z contains the transformation matrix produced by  eltran
  37.            after the reduction by  elmhes, or by  ortran  after the
  38.            reduction by  orthes, if performed.  if the eigenvectors
  39.            of the hessenberg matrix are desired, z must contain the
  40.            identity matrix.
  41.  
  42.       on output
  43.  
  44.          h has been destroyed.
  45.  
  46.          wr and wi contain the real and imaginary parts,
  47.            respectively, of the eigenvalues.  the eigenvalues
  48.            are unordered except that complex conjugate pairs
  49.            of values appear consecutively with the eigenvalue
  50.            having the positive imaginary part first.  if an
  51.            error exit is made, the eigenvalues should be correct
  52.            for indices ierr+1,...,n.
  53.  
  54.          z contains the real and imaginary parts of the eigenvectors.
  55.            if the i-th eigenvalue is real, the i-th column of z
  56.            contains its eigenvector.  if the i-th eigenvalue is complex
  57.            with positive imaginary part, the i-th and (i+1)-th
  58.            columns of z contain the real and imaginary parts of its
  59.            eigenvector.  the eigenvectors are unnormalized.  if an
  60.            error exit is made, none of the eigenvectors has been found.
  61.  
  62.          ierr is set to
  63.            zero       for normal return,
  64.            j          if the limit of 30*n iterations is exhausted
  65.                       while the j-th eigenvalue is being sought.
  66.  
  67.       calls cdiv for complex division.
  68.  
  69.  
  70.       this routine is a C-translation of the FORTRAN 77 source code
  71.       written by the mathematics and computer science division,
  72.       argonne national laboratory
  73.       last change :   september 1989.
  74.  
  75.       mark myers
  76.       Center for Applied Mathematics 
  77.       Cornell University    (607) 255-4195
  78.       ------------------------------------------------------------------ */
  79.  
  80.       ierr = 0;
  81.       norm = 0.e0;
  82.       k = 1;
  83.       for (i=1;i<=n;i++)                  /* store roots isolated bay balanc */
  84.         {for (j=k;j<=n;j++)               /* and compute matrix norm         */
  85.            norm = norm + fabs(h[i][j]);
  86.  
  87.          k = i;
  88.          if (i < *low | i>*igh)
  89.            {wr[i] = h[i][i];
  90.             wi[i] = 0.e0;} }
  91.  
  92.       en = *igh;    
  93.       t = 0.e0;
  94.       itn = 30*n;
  95.       while (en>= *low)           /* search for next eigenvalues */
  96.         {its = 0;
  97.          na = en - 1;
  98.          enm2 = na - 1;
  99. inner:   for (ll= *low;ll<=en;ll++) /* look for a single sub-diag element */
  100.            {l = en + *low - ll;   /* for l=en step -1 until low do... */
  101.             if (l == *low) 
  102.           break;
  103.             s = fabs(h[l-1][l-1]) + fabs(h[l][l]);
  104.             if (s == 0.e0)
  105.           s = norm;
  106.             tst1 = s;
  107.             tst2 = tst1 + fabs(h[l][l-1]);
  108.             if (tst2 == tst1)
  109.           break;}
  110.  
  111.          x = h[en][en];                      /* form shift */
  112.          if (l == en) 
  113.        {h[en][en] = x + t;       /* one root found */
  114.         wr[en] = h[en][en];
  115.         wi[en] = 0.e0;
  116.         en = na;
  117.         goto outer;}
  118.          y = h[na][na];
  119.          w = h[en][na] * h[na][en];
  120.          if (l == na) 
  121.        {p = (y - x) / 2.e0;      /* two roots found ... */
  122.             q = p * p + w;
  123.             zz = sqrt(fabs(q));
  124.             h[en][en] = x + t;
  125.             x = h[en][en];
  126.             h[na][na] = y + t;
  127.             if (q < 0.e0) 
  128.           {wr[na] = x + p;
  129.            wr[en] = x + p;
  130.            wi[na] = zz;
  131.            wi[en] = -zz;
  132.            en = enm2;
  133.            goto outer;}
  134.             zz = p + copysign(zz,p);  /* real roots */
  135.             wr[na] = x + zz;
  136.             wr[en] = wr[na];
  137.             if (zz != 0.e0) 
  138.            wr[en] = x - w / zz;
  139.             wi[na] = 0.e0;
  140.             wi[en] = 0.e0;
  141.             x = h[en][na];
  142.             s = fabs(x) + fabs(zz);
  143.             p = x / s;
  144.             q = zz / s;
  145.             r = sqrt(p*p+q*q);
  146.             p = p / r;
  147.             q = q / r;
  148.             for (j=na;j<=n;j++)       /* row modification */
  149.               {zz = h[na][j];
  150.                h[na][j] = q * zz + p * h[en][j];
  151.                h[en][j] = q * h[en][j] - p * zz;}
  152.             for (i=1;i<=en;i++)      /* column modification */
  153.               {zz = h[i][na];
  154.                h[i][na] = q * zz + p * h[i][en];
  155.                h[i][en] = q * h[i][en] - p * zz;}
  156.             for (i= *low;i<= *igh;i++)   /* accumulate transformations */ 
  157.               {zz = z[i][na];
  158.                z[i][na] = q * zz + p * z[i][en];
  159.                z[i][en] = q * z[i][en] - p * zz;}
  160.             en = enm2;
  161.         goto outer;}
  162.          if (itn == 0) 
  163.         return(en);
  164.          if (its == 10 | its == 20)
  165.            {t = t + x;    /* form exceptional shift */
  166.             for (i= *low;i<=en;i++)
  167.                h[i][i] = h[i][i] - x;
  168.             s = fabs(h[en][na]) + fabs(h[na][enm2]);
  169.             x = 0.75e0 * s;
  170.             y = x;
  171.             w = -0.4375e0 * s * s; }
  172.          its = its + 1;
  173.          itn = itn - 1; 
  174.  
  175.          for (mm=l;mm<=enm2;mm++)   /* look for 2 consecutive small */
  176.            {m = enm2 + l - mm;     /* sub-diagonal elements */
  177.             zz = h[m][m];          /* for m=en-2 step -1 until l do ... */
  178.             r = x - zz;
  179.             s = y - zz;
  180.             p = (r * s - w) / h[m+1][m] + h[m][m+1];
  181.             q = h[m+1][m+1] - zz - r - s;
  182.             r = h[m+2][m+1];
  183.             s = fabs(p) + fabs(q) + fabs(r);
  184.             p = p / s;
  185.             q = q / s;
  186.             r = r / s;
  187.             if (m == l)
  188.            break;
  189.             tst1 = fabs(p)*(fabs(h[m-1][m-1]) + fabs(zz) + fabs(h[m+1][m+1]));
  190.             tst2 = tst1 + fabs(h[m][m-1])*(fabs(q) + fabs(r));
  191.             if (tst2 == tst1) 
  192.            break;}
  193.  
  194.          mp2 = m + 2;
  195.  
  196.          for (i=mp2;i<=en;i++)
  197.            {h[i][i-2] = 0.e0;
  198.             if (i == mp2) 
  199.            continue;
  200.             h[i][i-3] = 0.e0;}
  201.          for (k=m;k<=na;k++)   /* double qr step involving rows l to en & */
  202.            {notlas = 'f';     /* columns m to en */
  203.         if (k != na)
  204.            notlas = 't';
  205.             if (k != m)
  206.           {p = h[k][k-1];
  207.                q = h[k+1][k-1];
  208.                r = 0.e0;
  209.                if (notlas == 't')
  210.               r = h[k+2][k-1];
  211.                x = fabs(p) + fabs(q) + fabs(r);
  212.                if (x == 0.e0) 
  213.               continue;  
  214.                p = p / x;
  215.                q = q / x;
  216.                r = r / x;}
  217.             s = copysign(sqrt(p*p+q*q+r*r),p);
  218.             if (k !=  m)
  219.                h[k][k-1] = -s * x; 
  220.             else
  221.                if (l != m)
  222.               h[k][k-1] = -h[k][k-1]; 
  223.             p = p + s;
  224.             x = p / s;
  225.             y = q / s;
  226.             zz = r / s;
  227.             q = q / p;
  228.             r = r / p;
  229.             if (notlas != 't')
  230.               {for (j=k;j<=n;j++)               /* row modification */
  231.              {p = h[k][j] + q * h[k+1][j];
  232.                   h[k][j] = h[k][j] - p * x;
  233.                   h[k+1][j] = h[k+1][j] - p * y;}
  234.                j = en;
  235.            if ( en>k+3 )
  236.               j = k+3;
  237.                for (i=1;i<=j;i++)              /* column modification */
  238.              {p = x * h[i][k] + y* h[i][k+1];
  239.               h[i][k] = h[i][k] - p;
  240.               h[i][k+1] = h[i][k+1] - p * q;}
  241.  
  242.                for (i= *low;i<= *igh;i++)  /* accumulate transformations */
  243.                  {p = x * z[i][k] + y * z[i][k+1];
  244.                   z[i][k] = z[i][k] - p;
  245.                   z[i][k+1] = z[i][k+1] - p * q;}    }
  246.             else
  247.               {for (j=k;j<=n;j++)        /* row modification */
  248.              {p = h[k][j] + q * h[k+1][j] + r * h[k+2][j];
  249.                   h[k][j] = h[k][j] - p * x;
  250.                   h[k+1][j] = h[k+1][j] - p * y;
  251.               h[k+2][j] = h[k+2][j] - p * zz;}
  252.                j = en;
  253.            if ( en>k+3 )
  254.               j = k+3;
  255.                for (i=1;i<=j;i++)       /* column modification */
  256.              {p = x * h[i][k] + y* h[i][k+1] + zz * h[i][k+2];
  257.               h[i][k] = h[i][k] - p;
  258.               h[i][k+1] = h[i][k+1] - p * q;
  259.               h[i][k+2] = h[i][k+2] - p * r;}
  260.  
  261.                for (i= *low;i<= *igh;i++)  /* accumulate transformations */
  262.                  {p = x * z[i][k] + y * z[i][k+1] + zz * z[i][k+2];
  263.                   z[i][k] = z[i][k] - p;
  264.                   z[i][k+1] = z[i][k+1] - p * q;
  265.                   z[i][k+2] = z[i][k+2] - p * r;} }
  266.            }
  267.              goto inner; 
  268.  
  269. outer:    ;}
  270.  
  271. /* all roots found.  backsubstitute to find vectors of upper triangular form */
  272.       if (norm == 0.e0)            
  273.     return;
  274.       for (nn=1;nn<=n;nn++)   /* for en = n step -1 until 1 do ..  */
  275.         {en = n + 1 - nn;
  276.          p = wr[en];
  277.          q = wi[en];
  278.          na = en - 1;
  279.      if (q > 0.e0)
  280.         goto dout;  
  281.          else if (q == 0.e0) 
  282.            {m = en;             /* real vector */ 
  283.             h[en][en] = 1.e0;
  284.             if (na == 0)
  285.            goto dout;
  286.  
  287.             for (ii=1;ii<=na;ii++)       /* for i=en-1 step -1 until 1 do .. */
  288.               {i = en - ii;
  289.                w = h[i][i] - p;
  290.                r = 0.e0;
  291.  
  292.                for (j=m;j<=en;j++)
  293.                   r = r + h[i][j] * h[j][en];
  294.  
  295.                if (wi[i] < 0.e0)
  296.                  {zz = w;
  297.                   s = r;
  298.           continue;}
  299.                m = i;
  300.                if (wi[i] == 0.e0)
  301.                  {t = w;
  302.                   if (t == 0.e0) 
  303.                     {tst1 = norm;
  304.                      t = tst1;
  305.                      while (tst2 > tst1)
  306.                        {t = 0.01e0 * t;
  307.                         tst2 = norm + t;}
  308.                     }
  309.                h[i][en] = -r / t;}
  310.             else
  311.               {x = h[i][i+1];       /* solve real equations */
  312.                y = h[i+1][i];
  313.                q = (wr[i] - p) * (wr[i] - p) + wi[i] * wi[i];
  314.                t = (x * s - zz * r) / q;
  315.                h[i][en] = t;
  316.                if (fabs(x) > fabs(zz)) 
  317.               h[i+1][en] = (-r - w * t) / x;
  318.                else
  319.                   h[i+1][en] = (-s - y * t) / zz;
  320.                t = fabs(h[i][en]);
  321.                if (t != 0.e0)
  322.                  {tst1 = t;
  323.                   tst2 = tst1 + 1.e0/tst1;
  324.                   if (tst2 <= tst1)
  325.                     for (j=i;j<=en;j++)
  326.                       h[j][en] = h[j][en]/t;}
  327.           }   }  }
  328.  
  329.        else if (q < 0.e0)        /* complex vector */
  330.      {m = na;
  331.  
  332. /* last vector component chosen imaginary so that
  333.    eigenvector matrix is triangular                */
  334.  
  335.          if (fabs(h[en][na]) > fabs(h[na][en]))
  336.            {h[na][na] = q / h[en][na];
  337.             h[na][en] = -(h[en][en] - p) / h[en][na];}
  338.          else       
  339.             cdiv(0.e0,-h[na][en],h[na][na]-p,q,&h[na][na],&h[na][en]);
  340.          h[en][na] = 0.e0;
  341.          h[en][en] = 1.e0;
  342.          enm2 = na - 1;
  343.          if (enm2 == 0) 
  344.         goto dout;
  345.          for (ii=1;ii<=enm2;ii++)     /* for i=en-2 step -1 until 1 do... */
  346.            {i = na - ii;
  347.             w = h[i][i] - p;
  348.             ra = 0.e0;
  349.             sa = 0.e0;
  350.             for (j=m;j<=en;j++) 
  351.               {ra = ra + h[i][j] * h[j][na];
  352.                sa = sa + h[i][j] * h[j][en];}
  353.  
  354.             if (wi[i] < 0.e0) 
  355.               {zz = w;
  356.                r = ra;
  357.                s = sa;
  358.                continue;}
  359.             m = i;
  360.             if (wi[i] == 0.e0)
  361.                cdiv(-ra,-sa,w,q,&h[i][na],&h[i][en]);
  362.             else      
  363.               {x = h[i][i+1];      /* solve complex equations */
  364.                y = h[i+1][i];
  365.                vr = (wr[i] - p) * (wr[i] - p) + wi[i] * wi[i] - q * q;
  366.                vi = (wr[i] - p) * 2.e0 * q;
  367.                if (vr == 0.e0 & vi == 0.e0)           
  368.                  {tst1 = norm * (fabs(w)+fabs(q)+fabs(x)+fabs(y)+fabs(zz));
  369.                   vr = tst1;
  370.           while (tst2 > tst1)
  371.                     {vr = 0.01e0 * vr;
  372.                      tst2 = tst1 + vr;}}
  373.                cdiv(x*r-zz*ra+q*sa,x*s-zz*sa-q*ra,vr,vi,&h[i][na],&h[i][en]);
  374.                if (fabs(x) > fabs(zz) + fabs(q)) 
  375.                 {h[i+1][na] = (-ra - w * h[i][na] + q * h[i][en]) / x;
  376.                  h[i+1][en] = (-sa - w * h[i][en] - q * h[i][na]) / x;}
  377.                else      
  378.                  cdiv(-r-y*h[i][na],-s-y*h[i][en],zz,q,&h[i+1][na],&h[i+1][en]);}
  379.             t = fabs(h[i][na]);
  380.         if (fabs(h[i][en]) > t)
  381.                t = fabs(h[i][en]);
  382.             if (t != 0.e0) 
  383.               {tst1 = t;
  384.                tst2 = tst1 + 1.e0/tst1;
  385.                if (tst2 <= tst1) 
  386.                  {for (j=i;j<=en;j++)
  387.                     {h[j][na] = h[j][na]/t;
  388.                      h[j][en] = h[j][en]/t;}}}
  389.                  
  390.         }   } 
  391.  
  392. dout:  ;
  393.      }
  394.   
  395.       for (i=1;i<=n;i++)
  396.         {if (i >= *low & i <= *igh) 
  397.        continue;
  398.          for (j=1;j<=n;j++)
  399.            z[i][j] = h[i][j];}
  400.  
  401.       for (jj= *low;jj<=n;jj++) /* multiply by transformation matrix to give */
  402.         {j = n + *low - jj;     /* vectors of original full matrix */
  403.          m = j;                /* for j=n step -1 until low do ... */
  404.          if (*igh < j)
  405.         m = *igh;
  406.          for (i= *low;i<= *igh;i++)
  407.            {zz = 0.e0;
  408.             for (k= *low;k<=m;k++)
  409.                zz = zz + z[i][k] * h[k][j];
  410.             z[i][j] = zz;}
  411.         }
  412. }
  413.